
#ifndef __GPU_2BWTSEARCH_CU__
#define __GPU_2BWTSEARCH_CU__

#include"2BWTSearch.h"
#include "Match.cu"
#include "mian_util.h"
#include "cuda_header.h"
#include "buffer.h"
#include <thrust/sort.h>

struct compuint3 {
	__host__ __device__ 
	bool operator()(uint3 a, uint3 b) {
		return a.x < b.x;
	}
};


void sortSA(uint3* cpuSA, uint32 numSA)
{
	thrust::sort(cpuSA, cpuSA+numSA, compuint3());
}

void BWT2GBWT(BWT *bwt, BWT **temp, unsigned int** d_occ)
{
	LookupTable * tempLookupTable = (LookupTable *)malloc(sizeof(LookupTable));

	BWT *d_bwt = (BWT *)malloc(sizeof(BWT));
	//	d_bwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	//	printf("size of BWT : %d\n", sizeof(BWT));

	unsigned int bwtCodeMemSize = bwt->bwtSizeInWord * sizeof(int);
	//     printf("bwtCodeMemSize : %d\n", bwtCodeMemSize);
	cutilSafeCall(cudaMalloc(&d_bwt->bwtCode, bwtCodeMemSize));
	cutilSafeCall(cudaMemcpy(d_bwt->bwtCode, bwt->bwtCode, bwtCodeMemSize, cudaMemcpyHostToDevice));

#if LAYER_INDEX_1
	unsigned int gpuOccValueMemSize = (ceil(bwt->textLength / (double)GPU_OCC_INTERVAL)* ALPHABET_SIZE * sizeof(int));
	cutilSafeCall(cudaMalloc(&d_bwt->GPUOccValue, gpuOccValueMemSize));
	cutilSafeCall(cudaMemcpy(d_bwt->GPUOccValue, bwt->GPUOccValue, gpuOccValueMemSize, cudaMemcpyHostToDevice));
#endif
#if LAYER_INDEX_2
	unsigned int bwtOccValueMemSize = bwt->occSizeInWord * sizeof(int);
	cutilSafeCall(cudaMalloc(&d_bwt->occValue,bwtOccValueMemSize));
	cutilSafeCall(cudaMemcpy(d_bwt->occValue, bwt->occValue, bwtOccValueMemSize, cudaMemcpyHostToDevice));

	unsigned int bwtOccMajorValueMemSize = bwt->occMajorSizeInWord * sizeof(int);
	cutilSafeCall(cudaMalloc(&d_bwt->occMajorValue, bwtOccMajorValueMemSize));
	cutilSafeCall(cudaMemcpy(d_bwt->occMajorValue, bwt->occMajorValue, bwtOccMajorValueMemSize, cudaMemcpyHostToDevice));

#endif
	if(NULL != d_occ) {
		(*d_occ) = d_bwt->GPUOccValue;
	}


	unsigned long long NR_TOP = 1 << ((bwt->lookupTable->tableSize) * 2);
	unsigned int lookupTableMemSize = NR_TOP * sizeof(int);
	//	     printf("lookupTableSize: %d\n", lookupTableMemSize);
	cutilSafeCall(cudaMalloc(&(tempLookupTable->table), lookupTableMemSize));
	cutilSafeCall(cudaMemcpy(tempLookupTable->table, bwt->lookupTable->table, lookupTableMemSize, cudaMemcpyHostToDevice));
	tempLookupTable->tableSize = bwt->lookupTable->tableSize;

	cutilSafeCall(cudaMalloc(&(d_bwt->lookupTable), sizeof(LookupTable)));
	cutilSafeCall(cudaMemcpy(d_bwt->lookupTable, tempLookupTable, sizeof(LookupTable), cudaMemcpyHostToDevice));


	d_bwt->saLength = bwt->saLength;
	d_bwt->textLength = bwt->textLength;
	//	d_bwt->lookupTableSize = bwt->lookupTable->Size;
	d_bwt->bwtSizeInWord = bwt->bwtSizeInWord;
	d_bwt->occSizeInWord = bwt->occSizeInWord;
	d_bwt->occMajorSizeInWord = bwt->occMajorSizeInWord;
	d_bwt->saValueSizeInWord = bwt->bwtSizeInWord;
	d_bwt->GPUOccValueInWord = ceil(bwt->textLength / (double)GPU_OCC_INTERVAL)* ALPHABET_SIZE;

	d_bwt->inverseSa0 = bwt->inverseSa0;
	d_bwt->saInterval = bwt->saInterval;


	unsigned int BWTMemSize = sizeof(BWT);
	//	printf("size of BWT : %d\n", BWTMemSize);
	cutilSafeCall(cudaMalloc(temp, BWTMemSize));
	cutilSafeCall(cudaMemcpy(*temp, d_bwt, BWTMemSize, cudaMemcpyHostToDevice));

	free(d_bwt);
	free(tempLookupTable);
}


void CudaFreeBWT(BWT *d_bwt)
{
	BWT *temp = (BWT *)malloc(sizeof(BWT));
	cutilSafeCall(cudaMemcpy(temp, d_bwt, sizeof(BWT), cudaMemcpyDeviceToHost));

	LookupTable * tempLookupTable =  (LookupTable *)malloc(sizeof(LookupTable));
	cutilSafeCall(cudaMemcpy(tempLookupTable, temp->lookupTable, sizeof(LookupTable), cudaMemcpyDeviceToHost));
	cutilSafeCall(cudaFree(tempLookupTable->table));

	cutilSafeCall(cudaFree(temp->bwtCode));
#if LAYER_INDEX_2
	cutilSafeCall(cudaFree(temp->occValue));
	cutilSafeCall(cudaFree(temp->occMajorValue));
#endif
#if LAYER_INDEX_1
	cutilSafeCall(cudaFree(temp->GPUOccValue));
#endif
	//	cutilSafeCall(cudaFree(temp->countValue));
	cutilSafeCall(cudaFree(d_bwt));

	free(temp);
	free(tempLookupTable);
}

void createOccValueFile(BWT *bwt, char* fileName)
{
	BWT * d_bwt;
	BWT2GBWT(bwt, &d_bwt);

	unsigned int occSize = ceil(bwt->textLength / (double)GPU_OCC_INTERVAL) * ALPHABET_SIZE;
	unsigned int * occValueFile = (unsigned int *) malloc(occSize * sizeof(int));
	unsigned int * d_occValueFile;

	cutilSafeCall(cudaMalloc(&d_occValueFile, occSize*sizeof(int)));
	cutilSafeCall(cudaMemcpy(d_occValueFile, occValueFile, occSize*sizeof(int), cudaMemcpyHostToDevice));

	CreateOccValue<<<BLOCK_SIZE, THREAD_SIZE>>>(d_bwt, d_occValueFile);

	cutilSafeCall(cudaMemcpy(occValueFile, d_occValueFile, occSize*sizeof(int), cudaMemcpyDeviceToHost));
	CudaFreeBWT(d_bwt);

	int i,j;
	for(i=0; i<100; i++){
		for(j=0; j<ALPHABET_SIZE; j++) {
			printf("%u\t", occValueFile[i*ALPHABET_SIZE + j]);
		}
		printf("\n");
	}

	FILE * fout = fopen(fileName,"wb");
	fwrite(occValueFile, sizeof(int), occSize, fout);
	fclose(fout);

}



void GPU2BWTInitIndex(BWT** d_bwt, BWT** d_revbwt, 
		BWT* bwt, BWT* revbwt, 
		st_MemoryBuf** GPUMemBuf, const uint64 gpuBufSize) {
	BWT2GBWT(bwt, d_bwt);
	BWT2GBWT(revbwt, d_revbwt);
	//*GPUMemBuf = new st_MemoryBuf(gpuBufSize, GMEM);

#if TEXTURE
	int texW = 65536;
	int texH = ceil(ceil(bwt->textLength / (double)GPU_OCC_INTERVAL) * ALPHABET_SIZE / (double)texW) ;  
	const int bytePitch = sizeof(int)*texW;
	BWT bwt_ptr , revbwt_ptr;
	cutilSafeCall(cudaMemcpy(&bwt_ptr, *d_bwt, sizeof(BWT), cudaMemcpyDeviceToHost));
	cutilSafeCall(cudaMemcpy(&revbwt_ptr, *d_revbwt, sizeof(BWT), cudaMemcpyDeviceToHost));

	cudaBindTexture2D(0, texOcc, bwt_ptr.GPUOccValue, texW, texH, bytePitch);
	cudaBindTexture2D(0, texRevOcc, revbwt_ptr.GPUOccValue, texW, texH, bytePitch);

	texOcc.normalized = 0;
	texRevOcc.normalized = 0;
#endif

	/*prepaer the constant memory*/
	unsigned int bwtCountValueMemSize = (ALPHABET_SIZE + 1) * sizeof(int);
	cutilSafeCall(cudaMemcpyToSymbol(CountValue, bwt->countValue, bwtCountValueMemSize,0, cudaMemcpyHostToDevice));
	int* h_shift = (int*)malloc(sizeof(int)*QUERY_LENGTH);
	for(int pos = 0; pos < QUERY_LENGTH; pos++) {
		h_shift[pos] = BIT_PER_CHAR*(CHAR_PER_WORD - 1  - (pos&(CHAR_PER_WORD - 1)));
	}
	cutilSafeCall(cudaMemcpyToSymbol(c_shiftTable, h_shift, sizeof(int)*QUERY_LENGTH));
	free(h_shift);	

	printf("GPU2BWTInit done.\n");
}




void GPU2BWTInit(BWT** d_bwt, BWT** d_revbwt, 
		BWT* bwt, BWT* revbwt, 
		st_MemoryBuf** GPUMemBuf, const uint64 gpuBufSize) {
	BWT2GBWT(bwt, d_bwt);
	BWT2GBWT(revbwt, d_revbwt);
	*GPUMemBuf = new st_MemoryBuf(gpuBufSize, GMEM);

#if TEXTURE
	int texW = 65536;
	int texH = ceil(ceil(bwt->textLength / (double)GPU_OCC_INTERVAL) * ALPHABET_SIZE / (double)texW) ;  
	const int bytePitch = sizeof(int)*texW;
	BWT bwt_ptr , revbwt_ptr;
	cutilSafeCall(cudaMemcpy(&bwt_ptr, *d_bwt, sizeof(BWT), cudaMemcpyDeviceToHost));
	cutilSafeCall(cudaMemcpy(&revbwt_ptr, *d_revbwt, sizeof(BWT), cudaMemcpyDeviceToHost));

	cudaBindTexture2D(0, texOcc, bwt_ptr.GPUOccValue, texW, texH, bytePitch);
	cudaBindTexture2D(0, texRevOcc, revbwt_ptr.GPUOccValue, texW, texH, bytePitch);

	texOcc.normalized = 0;
	texRevOcc.normalized = 0;
#endif

	/*prepaer the constant memory*/
	unsigned int bwtCountValueMemSize = (ALPHABET_SIZE + 1) * sizeof(int);
	cutilSafeCall(cudaMemcpyToSymbol(CountValue, bwt->countValue, bwtCountValueMemSize,0, cudaMemcpyHostToDevice));
	int* h_shift = (int*)malloc(sizeof(int)*QUERY_LENGTH);
	for(int pos = 0; pos < QUERY_LENGTH; pos++) {
		h_shift[pos] = BIT_PER_CHAR*(CHAR_PER_WORD - 1  - (pos&(CHAR_PER_WORD - 1)));
	}
	cutilSafeCall(cudaMemcpyToSymbol(c_shiftTable, h_shift, sizeof(int)*QUERY_LENGTH));
	free(h_shift);	

	printf("GPU2BWTInit done.\n");
}

uint64 GPU2BWTSearch1Kernel(const unsigned int * h_Queries, unsigned int queryNum,
		BWT *d_bwt, BWT *d_revbwt,
		st_MemoryBuf* GPUMemBuf) {

	PRINT_FUNC_NAME;

	unsigned int numQuery = queryNum;

	const uint64 gpuResultBufSize = 1700000000; //just used for result output per kernel
	uint3* d_result = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_result1 = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	int* d_numResultPerQuery = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQuery1 = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	unsigned int* d_queryBuf = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBuf1 = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* d_queryIdBuf = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);	
	int* d_queryIdBuf1 = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);	
	unsigned int* d_queryBuf2 = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* d_queryIdBuf2 = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);	
	int* d_temp = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);

#if SORT_QUERY  
	unsigned int * sortedQueryForward = SortQueriesForward(h_Queries, numQuery);
	unsigned int * sortedQueryBackward = SortQueriesBackward(h_Queries, numQuery);
	unsigned int* d_Queries_Fwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_Queries_Bwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries_Fwd, sortedQueryForward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries_Bwd, sortedQueryBackward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
#else
	unsigned int* d_Queries = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries, h_Queries, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
#endif
	int* origQueryId = (int*)malloc(sizeof(int)*queryNum);
	for(int i = 0; i < queryNum; i++) {
		origQueryId[i] = i;
	}
	int* d_QueryId = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	TOGPU(d_QueryId, origQueryId, sizeof(int)*queryNum);
	int* d_queryLabel = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	cutilSafeCall(cudaMemset(d_queryLabel, 0, sizeof(int)*queryNum));
	printf("GPU memory done.\n");

	int maxNumResultPerQuery;
	int numBlock, numThread;
	int sharedMemSize = 0;
	uint64 numAlignment = 0;



	printTitle("exact match");
	numBlock = 1204;
	numThread = 128;
#if SORT_QUERY
	exactMatch_kernel<<<numBlock, numThread>>>(d_result1,
			d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	exactMatch_kernel<<<numBlock, numThread>>>(d_result1,
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT exact match kernel 1");
	uint3* resultExact = (uint3 *)malloc(sizeof(uint3) * numQuery );
	cutilSafeCall(cudaMemcpy(resultExact, d_result1, sizeof(uint3)*numQuery, cudaMemcpyDeviceToHost));

	for(int i=0; i<numQuery; i++) {
		if(resultExact[i].x  < numQuery)
			numAlignment += resultExact[i].z - resultExact[i].y + 1;
	}
	printf("numAlignment exact match: %lu\n\n", numAlignment);


	/*******************************************************/
	/*1 mismatch*/
	/*1 mismatch forward*/
	//pass 1
	printTitle("1 mismatch");
	printf("1 mismatch forward search: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	numBlock = 1204;
	numThread = 128;
#if SORT_QUERY
	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries_Fwd, numQuery, d_bwt, d_revbwt);
#else
	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT1Mismatch forward_kernel 1");
	uint64 numAlignment1MisForwardPass1 = 0;
	int numUndone1MisForwardPass1 = 0;

#if SORT_QUERY
	collectUndoneGPU(       &numAlignment1MisForwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisForwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignment1MisForwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisForwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("1MisForward undone num = %lu\n", numUndone1MisForwardPass1);
	printf("numAlignment1MisForwardPass1 = %lu\n\n", numAlignment1MisForwardPass1);

	numAlignment += numAlignment1MisForwardPass1;

	//pass 2
	const int maxNumResultPerQuery1F = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisForwardPass1);

	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>
		(d_result1, d_numResultPerQuery1,
		 maxNumResultPerQuery1F, d_queryBuf1, numUndone1MisForwardPass1, d_bwt, d_revbwt);
	uint64 numAlignment1MisForwardPass2 = 0;
	int numUndone1MisForwardPass2 = 0;
	countAlignGPU(  &numAlignment1MisForwardPass2, d_result1, d_numResultPerQuery1,
			maxNumResultPerQuery1F, numUndone1MisForwardPass1,
			d_temp, d_queryIdBuf1, d_queryLabel);
	printf("numUndone1MisForwardPass2 = %d\n", numUndone1MisForwardPass2);
	assert(numUndone1MisForwardPass2 == 0);
	printf("numAlignment1MisForwardPass2 = %lu\n\n", numAlignment1MisForwardPass2);
	numAlignment += numAlignment1MisForwardPass2;



	/* backward */
	//pass 1
	printf("1 mismatch backward search: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	numBlock = 1204;
	numThread = 128;
#if SORT_QUERY
	BWT1MismatchBackward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries_Bwd, numQuery, d_bwt, d_revbwt);
#else
	BWT1MismatchBackward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT1Mismatch backward_kernel 1");


	uint64 numAlignment1MisBackwardPass1 = 0;
	int numUndone1MisBackwardPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignment1MisBackwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisBackwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Bwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignment1MisBackwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisBackwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("1MisBackward undone num = %lu\n", numUndone1MisBackwardPass1);
	printf("numAlignment1MisBackwardPass1 = %lu\n\n", numAlignment1MisBackwardPass1);
	numAlignment += numAlignment1MisBackwardPass1;
	//pass 2
	const int maxNumResultPerQuery1B = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisBackwardPass1);
	BWT1MismatchBackward_Kernel<<<numBlock, numThread, 0>>>
		(d_result1, d_numResultPerQuery1,
		 maxNumResultPerQuery1B, d_queryBuf1, numUndone1MisBackwardPass1, d_bwt, d_revbwt);
	uint64 numAlignment1MisBackwardPass2 = 0;
	int numUndone1MisBackwardPass2 = 0;
	countAlignGPU(  &numAlignment1MisBackwardPass2, d_result1, d_numResultPerQuery1,
			maxNumResultPerQuery1B, numUndone1MisBackwardPass1,
			d_temp, d_queryIdBuf1, d_queryLabel);
	printf("numUndone1MisBackwardPass2 = %d\n", numUndone1MisBackwardPass2);
	assert(numUndone1MisBackwardPass2 == 0);
	printf("numAlignment1MisBackwardPass2 = %lu\n\n", numAlignment1MisBackwardPass2);
	numAlignment += numAlignment1MisBackwardPass2;




	printTitle("2 mismatches");
	//case A, pass 1
	printf("4 Case kernel, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	printf("maxNumResultPerQuery: %u\n", maxNumResultPerQuery);	
	numBlock = 1024;
	numThread = 128;

#if SORT_QUERY
	BWT2MismatchSearch_4Case_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearch_4Case_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearch_4Case_kernel 1");
	uint64 numAlignmentPass1 = 0;
	int numUndonePass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignmentPass1, d_queryBuf, d_queryIdBuf, &numUndonePass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Bwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignmentPass1, d_queryBuf, d_queryIdBuf, &numUndonePass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignment 4Case kernel Pass1 = %lu\n\n", numAlignmentPass1);
	numAlignment += numAlignmentPass1;


	/********************* Pass2 *******************************************************/
	printf("4Case kernel, pass 2, #query = %d\n", numUndonePass1);
	uint64 numAlignmentPass2 = 0;
	int numUndonePass2 = 0;
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numUndonePass1);
	printf("maxNumResultPerQuery: %u\n", maxNumResultPerQuery);	

	BWT2MismatchSearch_4Case_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_queryBuf, d_queryIdBuf, numUndonePass1, d_bwt, d_revbwt);

	collectUndoneGPU(       &numAlignmentPass2, d_queryBuf2, d_queryIdBuf2, &numUndonePass2,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_queryBuf, d_queryIdBuf, numUndonePass1, d_temp, d_queryLabel);

	printf("numAlignment 4Case kernel Pass2 = %lu\n\n", numAlignmentPass2);
	numAlignment += numAlignmentPass2;

	/*********************    Pass 3 ***************************************************/

	printf("4Case kernel, pass 3, #query = %d\n", numUndonePass2);

	if(numUndonePass2 > 0) {
		uint64 numAlignmentPass3 = 0;
		maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numUndonePass2);
		printf("maxNumResultPerQuery: %u\n", maxNumResultPerQuery);

		BWT2MismatchSearch_4Case_kernel<<<numBlock, numThread, sharedMemSize>>>
			(d_result, d_numResultPerQuery,
			 maxNumResultPerQuery, d_queryBuf2, d_queryIdBuf2, numUndonePass2, d_bwt, d_revbwt);

		countAlignGPU(  &numAlignmentPass3, d_result, d_numResultPerQuery,
				maxNumResultPerQuery, numUndonePass2,
				d_temp, d_queryIdBuf2, d_queryLabel);
		printf("numAlignmentCaseAPass2 = %lu\n\n", numAlignmentPass3);
		numAlignment += numAlignmentPass3;	
	}
	/*
	   if(numUndonePass3 > 0) {
	   BWT2MismatchSearch_4Case_kernel<<<numBlock, numThread, sharedMemSize>>>
	   (d_result, d_numResultPerQuery,
	   maxNumResultPerQuery, d_queryBuf, d_queryIdBuf, numUndonePass2, d_bwt, d_revbwt);


	   }
	 */

	return numAlignment;
}



uint64 GPU2BWTSearch(const unsigned int * h_Queries, unsigned int queryNum, 
		BWT *d_bwt, BWT *d_revbwt,
		st_MemoryBuf* GPUMemBuf) {

	PRINT_FUNC_NAME;

	/*GPU memory*/
	unsigned int numQuery = queryNum;
	const uint64 gpuResultBufSize = 500000000; //just used for result output per kernel
	uint3* d_result = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	int* d_numResultPerQuery = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	///now for the concurrent kernels
	unsigned int* d_queryBufA = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBufB = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBufC = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBufD = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBuf1 = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* d_queryIdBufA = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBufB = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBufC = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBufD = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBuf1 = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);


	uint3* d_resultA = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_resultB = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_resultC = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_resultD = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_result1 = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	int* d_numResultPerQueryA = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQueryB = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQueryC = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQueryD = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQuery1 = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_temp = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);


	/***********************Sort query***********************/
#if SORT_QUERY	
	unsigned int * sortedQueryForward = SortQueriesForward(h_Queries, numQuery);
	unsigned int * sortedQueryBackward = SortQueriesBackward(h_Queries, numQuery);
	unsigned int* d_Queries_Fwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_Queries_Bwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries_Fwd, sortedQueryForward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries_Bwd, sortedQueryBackward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
#else	
	unsigned int* d_Queries = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries, h_Queries, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
#endif
	int* origQueryId = (int*)malloc(sizeof(int)*queryNum);
	for(int i = 0; i < queryNum; i++) {
		origQueryId[i] = i;
	}
	int* d_QueryId = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	TOGPU(d_QueryId, origQueryId, sizeof(int)*queryNum);
	int* d_queryLabel = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	cutilSafeCall(cudaMemset(d_queryLabel, 0, sizeof(int)*queryNum));
	printf("GPU memory done.\n");



	/*kernel configuration and others*/
#ifdef PREFER_L1
	cudaFuncSetCacheConfig(exactMatch_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferL1);
#else 
	cudaFuncSetCacheConfig(exactMatch_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferShared);

#endif
	int maxNumResultPerQuery;
	int numBlock, numThread;
	int sharedMemSize = 0;
	uint64 numAlignment = 0;
	int nstreams = 6;
	cudaStream_t *streams = (cudaStream_t*) malloc(nstreams * sizeof(cudaStream_t));
	for(int i = 0; i < nstreams; i++) {
		cutilSafeCall( cudaStreamCreate(&(streams[i])) );
	}
	numBlock = 3584;
	numThread = 64;
	sharedMemSize = 0; //sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numThread;


	/*******************************************************/
	/*exact match pass1 only*/
	printTitle("exact match");
	numBlock = 1204;
	numThread = 128;
#if SORT_QUERY
	exactMatch_kernel<<<numBlock, numThread>>>(d_result1, 
			d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	exactMatch_kernel<<<numBlock, numThread>>>(d_result1, 
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT exact match kernel 1");
	uint3* resultExact = (uint3 *)malloc(sizeof(uint3) * numQuery );
	cutilSafeCall(cudaMemcpy(resultExact, d_result1, sizeof(uint3)*numQuery, cudaMemcpyDeviceToHost));

	for(int i=0; i<numQuery; i++) {
		if(resultExact[i].x  < numQuery)
			numAlignment += resultExact[i].z - resultExact[i].y + 1;	
	}
	printf("numAlignment exact match: %lu\n\n", numAlignment);


	/*******************************************************/
	/*1 mismatch*/
	/*1 mismatch forward*/
	//pass 1
	printTitle("1 mismatch");
	printf("1 mismatch forward search: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	numBlock = 1204;
	numThread = 128;
#if SORT_QUERY
	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries_Fwd, numQuery, d_bwt, d_revbwt);
#else
	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries, numQuery, d_bwt, d_revbwt);
#endif	
	cutilCheckMsg("BWT1Mismatch forward_kernel 1");
	uint64 numAlignment1MisForwardPass1 = 0;
	int numUndone1MisForwardPass1 = 0;

#if SORT_QUERY
	collectUndoneGPU(       &numAlignment1MisForwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisForwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else	
	collectUndoneGPU(       &numAlignment1MisForwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisForwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif	
	printf("1MisForward undone num = %lu\n", numUndone1MisForwardPass1);
	printf("numAlignment1MisForwardPass1 = %lu\n\n", numAlignment1MisForwardPass1);

	numAlignment += numAlignment1MisForwardPass1;
	//pass 2
	const int maxNumResultPerQuery1F = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisForwardPass1);

	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>
		(d_result1, d_numResultPerQuery1,
		 maxNumResultPerQuery1F, d_queryBuf1, numUndone1MisForwardPass1, d_bwt, d_revbwt);
	uint64 numAlignment1MisForwardPass2 = 0;
	int numUndone1MisForwardPass2 = 0;
	countAlignGPU(  &numAlignment1MisForwardPass2, d_result1, d_numResultPerQuery1,
			maxNumResultPerQuery1F, numUndone1MisForwardPass1,
			d_temp, d_queryIdBuf1, d_queryLabel);
	printf("numUndone1MisForwardPass2 = %d\n", numUndone1MisForwardPass2);
	assert(numUndone1MisForwardPass2 == 0);
	printf("numAlignment1MisForwardPass2 = %lu\n\n", numAlignment1MisForwardPass2);
	numAlignment += numAlignment1MisForwardPass2;



	/* backward */
	//pass 1
	printf("1 mismatch backward search: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	numBlock = 1204;
	numThread = 128;
#if SORT_QUERY
	BWT1MismatchBackward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries_Bwd, numQuery, d_bwt, d_revbwt);
#else
	BWT1MismatchBackward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT1Mismatch backward_kernel 1");


	uint64 numAlignment1MisBackwardPass1 = 0;
	int numUndone1MisBackwardPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignment1MisBackwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisBackwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Bwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignment1MisBackwardPass1, d_queryBuf1, d_queryIdBuf1, &numUndone1MisBackwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif	
	printf("1MisBackward undone num = %lu\n", numUndone1MisBackwardPass1);
	printf("numAlignment1MisBackwardPass1 = %lu\n\n", numAlignment1MisBackwardPass1);
	numAlignment += numAlignment1MisBackwardPass1;
	//pass 2
	const int maxNumResultPerQuery1B = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisBackwardPass1);
	BWT1MismatchBackward_Kernel<<<numBlock, numThread, 0>>>
		(d_result1, d_numResultPerQuery1,
		 maxNumResultPerQuery1B, d_queryBuf1, numUndone1MisBackwardPass1, d_bwt, d_revbwt);
	uint64 numAlignment1MisBackwardPass2 = 0;
	int numUndone1MisBackwardPass2 = 0;
	countAlignGPU(  &numAlignment1MisBackwardPass2, d_result1, d_numResultPerQuery1,
			maxNumResultPerQuery1B, numUndone1MisBackwardPass1,
			d_temp, d_queryIdBuf1, d_queryLabel);
	printf("numUndone1MisBackwardPass2 = %d\n", numUndone1MisBackwardPass2);
	assert(numUndone1MisBackwardPass2 == 0);
	printf("numAlignment1MisBackwardPass2 = %lu\n\n", numAlignment1MisBackwardPass2);
	numAlignment += numAlignment1MisBackwardPass2;





	/*******************************************************/
	/*two mismatches*/
	printTitle("2 mismatches");
	//case A, pass 1
	printf("Case A, pass 1, numQuery: %d\n", numQuery);
#if OPT_NUM_RESULT
	maxNumResultPerQuery = OPT_NUM_RESULT;
#else
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#endif
#if SORT_QUERY
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseA_kernel 1");
	uint64 numAlignmentCaseAPass1 = 0;
	int numUndoneCaseAPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(	&numAlignmentCaseAPass1, d_queryBufA, d_queryIdBufA, &numUndoneCaseAPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries_Bwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(	&numAlignmentCaseAPass1, d_queryBufA, d_queryIdBufA, &numUndoneCaseAPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseAPass1 = %lu\n\n", numAlignmentCaseAPass1);
	numAlignment += numAlignmentCaseAPass1;



	//case B, pass 1
	numQuery = queryNum;
	printf("Case B, pass 1, numQuery: %d\n", numQuery);
#if OPT_NUM_RESULT
	maxNumResultPerQuery = OPT_NUM_RESULT;
#else
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#endif

#if SORT_QUERY
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries_Fwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseB_kernel 1");
	uint64 numAlignmentCaseBPass1 = 0;
	int numUndoneCaseBPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(	&numAlignmentCaseBPass1, d_queryBufB, d_queryIdBufB, &numUndoneCaseBPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(	&numAlignmentCaseBPass1, d_queryBufB, d_queryIdBufB, &numUndoneCaseBPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseBPass1 = %lu\n\n", numAlignmentCaseBPass1);
	numAlignment += numAlignmentCaseBPass1;


	//case C, pass 1
	numQuery = queryNum;
	printf("Case C, pass 1, numQuery: %d\n", numQuery);
#if OPT_NUM_RESULT
	maxNumResultPerQuery = OPT_NUM_RESULT;
#else
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#endif

#if SORT_QUERY
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries_Fwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseC_kernel 1");
	uint64 numAlignmentCaseCPass1 = 0;
	int numUndoneCaseCPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(	&numAlignmentCaseCPass1, d_queryBufC, d_queryIdBufC, &numUndoneCaseCPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(	&numAlignmentCaseCPass1, d_queryBufC, d_queryIdBufC, &numUndoneCaseCPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseCPass1 = %lu\n\n", numAlignmentCaseCPass1);
	numAlignment += numAlignmentCaseCPass1;


	//case D, pass 1
	numQuery = queryNum;
	printf("Case D, pass 1, numQuery: %d\n", numQuery);
#if OPT_NUM_RESULT
	maxNumResultPerQuery = OPT_NUM_RESULT;
#else
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#endif

#if SORT_QUERY
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries_Fwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseD_kernel 1");
	uint64 numAlignmentCaseDPass1 = 0;
	int numUndoneCaseDPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(	&numAlignmentCaseDPass1, d_queryBufD, d_queryIdBufD, &numUndoneCaseDPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(	&numAlignmentCaseDPass1, d_queryBufD, d_queryIdBufD, &numUndoneCaseDPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseDPass1 = %lu\n\n", numAlignmentCaseDPass1);
	numAlignment += numAlignmentCaseDPass1;



	/////////////////////////////////////////////////////
	//Pass 2
	numBlock = 448;
	numThread = 64;
	sharedMemSize = 0; //sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numThread;

	/////////////////////
	//pass 2, stream execution
	printf("pass 2: numBlock = %d, numThread = %d\n", numBlock, numThread);
	const int maxNumResultPerQueryA2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseAPass1);
	const int maxNumResultPerQueryB2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseBPass1);
	const int maxNumResultPerQueryC2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseCPass1);
	const int maxNumResultPerQueryD2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseDPass1);





#ifdef CONCURRENT_EXEC
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize, streams[0]>>>
		(d_resultA, d_numResultPerQueryA, 
		 maxNumResultPerQueryA2, d_queryBufA, d_queryIdBufA, 
		 numUndoneCaseAPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize, streams[1]>>>
		(d_resultB, d_numResultPerQueryB, 
		 maxNumResultPerQueryB2, d_queryBufB, d_queryIdBufB, 
		 numUndoneCaseBPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize, streams[2]>>>
		(d_resultC, d_numResultPerQueryC, 
		 maxNumResultPerQueryC2, d_queryBufC, d_queryIdBufC, 
		 numUndoneCaseCPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize, streams[3]>>>
		(d_resultD, d_numResultPerQueryD, 
		 maxNumResultPerQueryD2, d_queryBufD, d_queryIdBufD, 
		 numUndoneCaseDPass1, d_bwt, d_revbwt);
#else 
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_resultA, d_numResultPerQueryA,
		 maxNumResultPerQueryA2, d_queryBufA, d_queryIdBufA,
		 numUndoneCaseAPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_resultB, d_numResultPerQueryB,
		 maxNumResultPerQueryB2, d_queryBufB, d_queryIdBufB,
		 numUndoneCaseBPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_resultC, d_numResultPerQueryC,
		 maxNumResultPerQueryC2, d_queryBufC, d_queryIdBufC,
		 numUndoneCaseCPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_resultD, d_numResultPerQueryD,
		 maxNumResultPerQueryD2, d_queryBufD, d_queryIdBufD,
		 numUndoneCaseDPass1, d_bwt, d_revbwt);
#endif

	///////////////////////////////////////
	printf("Case A, pass 2, #query = %d\n", numUndoneCaseAPass1);
	uint64 numAlignmentCaseAPass2 = 0;
	countAlignGPU(	&numAlignmentCaseAPass2, d_resultA, d_numResultPerQueryA, 
			maxNumResultPerQueryA2, numUndoneCaseAPass1,
			d_temp, d_queryIdBufA, d_queryLabel);
	printf("numAlignmentCaseAPass2 = %lu\n\n", numAlignmentCaseAPass2);
	numAlignment += numAlignmentCaseAPass2;



	printf("Case B, pass 2, #query = %d\n", numUndoneCaseBPass1);
	uint64 numAlignmentCaseBPass2 = 0;
	countAlignGPU(	&numAlignmentCaseBPass2, d_resultB, d_numResultPerQueryB, 
			maxNumResultPerQueryB2, numUndoneCaseBPass1,
			d_temp, d_queryIdBufB, d_queryLabel);
	printf("numAlignmentCaseBPass2 = %lu\n\n", numAlignmentCaseBPass2);
	numAlignment += numAlignmentCaseBPass2;



	printf("Case C, pass 2, #query = %d\n", numUndoneCaseCPass1);
	uint64 numAlignmentCaseCPass2 = 0;
	countAlignGPU(	&numAlignmentCaseCPass2, d_resultC, d_numResultPerQueryC, 
			maxNumResultPerQueryC2, numUndoneCaseCPass1,
			d_temp, d_queryIdBufC, d_queryLabel);
	printf("numAlignmentCaseCPass2 = %lu\n\n", numAlignmentCaseCPass2);
	numAlignment += numAlignmentCaseCPass2;




	printf("Case D, pass 2, #query = %d\n", numUndoneCaseDPass1);
	uint64 numAlignmentCaseDPass2 = 0;
	countAlignGPU(	&numAlignmentCaseDPass2, d_resultD, d_numResultPerQueryD, 
			maxNumResultPerQueryD2, numUndoneCaseDPass1,
			d_temp, d_queryIdBufD, d_queryLabel);
	printf("numAlignmentCaseDPass2 = %lu\n\n", numAlignmentCaseDPass2);
	numAlignment += numAlignmentCaseDPass2;


	const int numMatched = getNumMatchedQuery(d_queryLabel, queryNum);
	printf("===> numAlignment: %lu\n", numAlignment);
	printf("===> #matched query: %d\n", numMatched);

#if SORT_QUERY
	free(sortedQueryForward);
	free(sortedQueryBackward);
#endif
	return numAlignment;
}

uint64 GPU2BWTSearch_shareBuf(const unsigned int * h_Queries, unsigned int queryNum,
		BWT *d_bwt, BWT *d_revbwt,
		st_MemoryBuf* GPUMemBuf, uint3* h_result, uint32* h_resultNum)
{
	PRINT_FUNC_NAME; 
	/*GPU memory*/
	unsigned int numQuery = queryNum;
	const uint64 gpuResultBufSize = 2400000000; //just used for result output per kernel 
	uint3* d_result = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	int2* d_numResultPerQuery = (int2*)GPUMemBuf->newMem(sizeof(int2)* numQuery);
	unsigned int* d_queryBuf = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* d_queryIdBuf = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	unsigned int* d_queryBuf_temp = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* d_queryIdBuf_temp = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_temp = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	
	unsigned int *d_queryBufExec;
	int *d_queryIdBufExec;

	/***********************Sort query***********************/
#if SORT_QUERY   
	unsigned int * sortedQueryForward = SortQueriesForward(h_Queries, numQuery);
	unsigned int * sortedQueryBackward = SortQueriesBackward(h_Queries, numQuery);
	unsigned int* d_Queries_Fwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_Queries_Bwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries_Fwd, sortedQueryForward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries_Bwd, sortedQueryBackward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
#else
	unsigned int* d_Queries = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries, h_Queries, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
#endif
	int* origQueryId = (int*)malloc(sizeof(int)*queryNum);
	for(int i = 0; i < queryNum; i++) {
		origQueryId[i] = i;
	}
	int* d_QueryId = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	TOGPU(d_QueryId, origQueryId, sizeof(int)*queryNum);
	int* d_queryLabel = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	cutilSafeCall(cudaMemset(d_queryLabel, 0, sizeof(int)*queryNum));
	printf("GPU memory done.\n");

	/*kernel configuration and others*/
#ifdef PREFER_L1 
	cudaFuncSetCacheConfig(exactMatch_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferL1);
#else
	cudaFuncSetCacheConfig(exactMatch_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferShared);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferShared);

#endif
	int maxNumResultPerQuery;
	int maxNumResultPerThread; 		
	int numBlock, numThread, passId;
	int sharedMemSize = 0;
	uint64 numAlignment = 0;
	uint64 numUndone = 0;
	uint32 numResult = 0;
	sharedMemSize = 0; //sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numThread; 

	/*******************************************************/ 
	/*exact match pass1 only*/ 
	printTitle("exact match"); 
	numBlock = 512; 
	//numBlock = 1204; 
	numThread = 128; 
#if SORT_QUERY 
	exactMatch_kernel<<<numBlock, numThread>>>(d_result,  
			d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt); 
#else 
	exactMatch_kernel<<<numBlock, numThread>>>(d_result,  
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt); 
#endif 
	cutilCheckMsg("BWT exact match kernel 1"); 
//	uint3* resultExact = (uint3 *)malloc(sizeof(uint3) * numQuery ); 
//	cutilSafeCall(cudaMemcpy(resultExact, d_result, sizeof(uint3)*numQuery, cudaMemcpyDeviceToHost)); 


/*	int count = 0;
	for(int i=0; i<numQuery; i++) { 
		if(resultExact[i].x  < numQuery) {
			numAlignment += resultExact[i].z - resultExact[i].y + 1;         
			count++;
		}
	} 
	printf("count = %d\n", count);

	printf("numAlignment exact match: %lu\n\n", numAlignment); 
*/


	collectExactResult(d_result, numQuery, h_result, h_resultNum);	
	printf("*h_resultNum = %d\n", *h_resultNum);

	for(int i=0; i<(*h_resultNum); i++)
		numAlignment += h_result[i].z - h_result[i].y + 1;
		
	printf("numAlignment exact match: %u\n\n", (numAlignment)); 


	/*******************************************************/ 
	/*1 mismatch*/ 
	/*1 mismatch forward*/ 
	//pass 1 
	printTitle("1 mismatch"); 
	printf("1 mismatch forward search: %d\n", numQuery); 
	numBlock = 512; 
	numThread = 128; 
	passId = 0;
	numUndone = numQuery;
#if SORT_QUERY
	d_queryBufExec = d_Queries_Fwd;
#else
	d_queryBufExec = d_Queries;
#endif
	d_queryIdBufExec = d_QueryId;
	while(numUndone > 0) {
		if(numUndone < SMALL_DATA_SIZE){
			numBlock = SMALL_BLOCK_SIZE;
			numThread = SMALL_THREAD_SIZE;
		}
		maxNumResultPerThread = getNumOutputPerThread(gpuResultBufSize, numBlock, numThread, numUndone);
		cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
		cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
		passId++;
		sharedMemSize = numThread * OCC_COUNT_PER_THREAD * sizeof(int);
		BWT1MismatchForward_shareBuf_kernel<<<numBlock, numThread, sharedMemSize>>>(d_result, 
				d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, numUndone, d_bwt, d_revbwt); 
		cutilCheckMsg("BWT1Mismatch forward_kernel 1"); 
		uint64 numAlignment1MisForwardPass1 = 0; 
		int numUndone1MisForwardPass1 = 0; 

		collectUndoneGPU_shareBuf(       &numAlignment1MisForwardPass1, 
				d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
				d_result, d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_temp, d_queryLabel, numBlock,numThread, h_result, h_resultNum); 
	
		printf("pass %d 1MisForward undone num = %lu\n",passId, numUndone1MisForwardPass1); 
		printf("numAlignment1MisForward = %lu\n\n", numAlignment1MisForwardPass1); 

		numAlignment += numAlignment1MisForwardPass1; 
		numUndone = numUndone1MisForwardPass1;
		
		cutilSafeCall(cudaMemcpy(d_queryBuf_temp,d_queryBuf, 
				numQuery * QUERY_SIZE_IN_WORD * sizeof(int),cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemcpy(d_queryIdBuf_temp,d_queryIdBuf,
				numQuery*sizeof(int), cudaMemcpyDeviceToDevice));
	
		cutilSafeCall(cudaMemset(d_queryBuf, 0, numQuery*QUERY_SIZE_IN_WORD*sizeof(int)));
		cutilSafeCall(cudaMemset(d_queryIdBuf, 0, numQuery*sizeof(int)));
	
		d_queryBufExec = d_queryBuf_temp;
		d_queryIdBufExec = d_queryIdBuf_temp;
	}
	/*
	//pass 2 
	if(numUndone1MisForwardPass1 > 0) {
	numBlock = 128;
	numThread = 32;
	cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
	cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisForwardPass1); 
	maxNumResultPerThread = maxNumResultPerQuery * floor(numUndone1MisForwardPass1 / (double)(numBlock * numThread));
	if(maxNumResultPerThread == 0)
	maxNumResultPerThread = maxNumResultPerQuery;

	BWT1MismatchForward_shareBuf_kernel<<<numBlock, numThread>>> 
	(d_result, d_numResultPerQuery, 
	maxNumResultPerQuery, d_queryBuf, numUndone1MisForwardPass1, d_bwt, d_revbwt); 

	uint64 numAlignment1MisForwardPass2 = 0; 
	int numUndone1MisForwardPass2 = 0; 
	countAlignGPU_shareBuf(  &numAlignment1MisForwardPass2, d_result, d_numResultPerQuery, 
	maxNumResultPerThread, numUndone1MisForwardPass1, 
	d_temp, d_queryIdBuf, d_queryLabel,numBlock, numThread); 
	printf("numUndone1MisForwardPass2 = %d\n", numUndone1MisForwardPass2); 
	assert(numUndone1MisForwardPass2 == 0); 
	printf("numAlignment1MisForwardPass2 = %lu\n\n", numAlignment1MisForwardPass2); 
	numAlignment += numAlignment1MisForwardPass2; 
	}
	 */
	/* backward */ 
	//pass 1 
	printf("1 mismatch backward search: %d\n", numQuery);

	numBlock = 512; 
	//      numBlock = 1204; 
	numThread = 128;
	passId = 0;
	numUndone = numQuery;
#if SORT_QUERY
	d_queryBufExec = d_Queries_Bwd;
#else
	d_queryBufExec = d_Queries;
#endif
	d_queryIdBufExec = d_QueryId;


	while(numUndone > 0) {	
	                if(numUndone < SMALL_DATA_SIZE){
                        numBlock = SMALL_BLOCK_SIZE;
                        numThread = SMALL_THREAD_SIZE;
                }
		maxNumResultPerThread = getNumOutputPerThread(gpuResultBufSize, numBlock, numThread, numQuery);
		cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
		cutilSafeCall(cudaMemset(d_result, 0,gpuResultBufSize));
		passId++;
		sharedMemSize = numThread * OCC_COUNT_PER_THREAD * sizeof(int);
		BWT1MismatchBackward_shareBuf_kernel<<<numBlock, numThread, sharedMemSize>>>(d_result, 
				d_numResultPerQuery,maxNumResultPerThread, 
				d_queryBufExec, numUndone, d_bwt, d_revbwt); 

		cutilCheckMsg("BWT1Mismatch backward_kernel 1"); 

		uint64 numAlignment1MisBackwardPass1 = 0; 
		int numUndone1MisBackwardPass1 = 0; 
		collectUndoneGPU_shareBuf(&numAlignment1MisBackwardPass1, d_queryBuf, 
				d_queryIdBuf, &numUndone1MisBackwardPass1, 
				d_result, d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_temp, d_queryLabel, numBlock, numThread, h_result, h_resultNum); 
		printf("pass %d 1MisBackward undone num = %lu\n",passId, numUndone1MisBackwardPass1); 
		printf("numAlignment1MisBackward = %lu\n\n", numAlignment1MisBackwardPass1); 
		numAlignment += numAlignment1MisBackwardPass1; 
		
		numUndone = numUndone1MisBackwardPass1;
		cutilSafeCall(cudaMemcpy(d_queryBuf_temp,d_queryBuf, 
				numQuery * QUERY_SIZE_IN_WORD * sizeof(int),cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemcpy(d_queryIdBuf_temp,d_queryIdBuf,
				numQuery*sizeof(int), cudaMemcpyDeviceToDevice));
		
		cutilSafeCall(cudaMemset(d_queryBuf, 0, numQuery*QUERY_SIZE_IN_WORD*sizeof(int)));
		cutilSafeCall(cudaMemset(d_queryIdBuf, 0, numQuery*sizeof(int)));
		d_queryBufExec = d_queryBuf_temp;
		d_queryIdBufExec = d_queryIdBuf_temp;
	}
	/*
	//pass 2 
	if(numUndone1MisBackwardPass1 > 0) {
	numBlock = 128;
	numThread = 32;
	cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
	cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
	const int maxNumResultPerQuery1B = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisBackwardPass1); 
	maxNumResultPerThread = maxNumResultPerQuery1B * floor(numUndone1MisBackwardPass1 / (double)(numBlock * numThread));
	if(maxNumResultPerThread == 0)
	maxNumResultPerThread = maxNumResultPerQuery1B;
	BWT1MismatchBackward_shareBuf_kernel<<<numBlock, numThread, 0>>> 
	(d_result, d_numResultPerQuery, 
	maxNumResultPerQuery1B, d_queryBuf, numUndone1MisBackwardPass1, d_bwt, d_revbwt); 
	uint64 numAlignment1MisBackwardPass2 = 0; 
	int numUndone1MisBackwardPass2 = 0; 
	countAlignGPU_shareBuf(  &numAlignment1MisBackwardPass2, d_result, d_numResultPerQuery, 
	maxNumResultPerThread, numUndone1MisBackwardPass1, 
	d_temp, d_queryIdBuf, d_queryLabel, numBlock, numThread); 
	printf("numUndone1MisBackwardPass2 = %d\n", numUndone1MisBackwardPass2); 
	assert(numUndone1MisBackwardPass2 == 0); 
	printf("numAlignment1MisBackwardPass2 = %lu\n\n", numAlignment1MisBackwardPass2); 
	numAlignment += numAlignment1MisBackwardPass2; 
	} 
	 */

	/*******************************************************/ 
	/*two mismatches*/ 
	printTitle("2 mismatches"); 
	numBlock = 512; 
	numThread = 128; 
	passId = 0;
	numUndone = numQuery;
#if SORT_QUERY
	d_queryBufExec = d_Queries_Bwd;
#else
	d_queryBufExec = d_Queries;
#endif
	d_queryIdBufExec = d_QueryId;
	while(numUndone > 0) {
		if(numUndone < SMALL_DATA_SIZE){
			numBlock = SMALL_BLOCK_SIZE;
			numThread = SMALL_THREAD_SIZE;
		}
		maxNumResultPerThread = getNumOutputPerThread(gpuResultBufSize, numBlock, numThread, numUndone);
		cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
		cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
		passId++;
		BWT2MismatchSearchCaseA_shareBuf_kernel<<<numBlock, numThread>>>(d_result, 
				d_numResultPerQuery,maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_bwt, d_revbwt); 
		cutilCheckMsg("BWT1Mismatch forward_kernel 1"); 
		
		uint64 numAlignment1MisForwardPass1 = 0; 
		int numUndone1MisForwardPass1 = 0; 

		collectUndoneGPU_shareBuf(       &numAlignment1MisForwardPass1, 
				d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
				d_result, d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_temp, d_queryLabel, numBlock,numThread, h_result, h_resultNum); 
		
		printf("pass %d caseA undone num = %lu\n",passId, numUndone1MisForwardPass1); 
		printf("numAlignmentCaseA= %lu\n\n", numAlignment1MisForwardPass1); 

		numAlignment += numAlignment1MisForwardPass1; 
		numUndone = numUndone1MisForwardPass1;
		
		cutilSafeCall(cudaMemcpy(d_queryBuf_temp,d_queryBuf, 
				numQuery * QUERY_SIZE_IN_WORD * sizeof(int),cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemcpy(d_queryIdBuf_temp,d_queryIdBuf,
				numQuery*sizeof(int), cudaMemcpyDeviceToDevice));

		cutilSafeCall(cudaMemset(d_queryBuf, 0, numQuery*QUERY_SIZE_IN_WORD*sizeof(int)));
		cutilSafeCall(cudaMemset(d_queryIdBuf, 0, numQuery*sizeof(int)));
                d_queryBufExec = d_queryBuf_temp;
                d_queryIdBufExec = d_queryIdBuf_temp;
	}

	numBlock = 512; 
	numThread = 128; 
	passId = 0;
	numUndone = numQuery;
#if SORT_QUERY
	d_queryBufExec = d_Queries_Fwd;
#else
	d_queryBufExec = d_Queries;
#endif
	d_queryIdBufExec = d_QueryId;
	while(numUndone > 0) {
		if(numUndone < SMALL_DATA_SIZE){
			numBlock = SMALL_BLOCK_SIZE;
			numThread = SMALL_THREAD_SIZE;
		}
		maxNumResultPerThread = getNumOutputPerThread(gpuResultBufSize, numBlock, numThread, numUndone);
		cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
		cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
		passId++;
		sharedMemSize = numThread * OCC_COUNT_PER_THREAD * sizeof(int);
		BWT2MismatchSearchCaseB_shareBuf_kernel<<<numBlock, numThread, sharedMemSize>>>(d_result, 
				d_numResultPerQuery,maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_bwt, d_revbwt); 
		cutilCheckMsg("BWT1Mismatch forward_kernel 1"); 
		
		uint64 numAlignment1MisForwardPass1 = 0; 
		int numUndone1MisForwardPass1 = 0; 

		collectUndoneGPU_shareBuf(       &numAlignment1MisForwardPass1, 
				d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
				d_result, d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_temp, d_queryLabel, numBlock,numThread, h_result, h_resultNum); 
		
		printf("pass %d caseB undone num = %lu\n",passId, numUndone1MisForwardPass1); 
		printf("numAlignmentCaseB = %lu\n\n", numAlignment1MisForwardPass1); 

		numAlignment += numAlignment1MisForwardPass1; 
		numUndone = numUndone1MisForwardPass1;
		cutilSafeCall(cudaMemcpy(d_queryBuf_temp,d_queryBuf, 
				numQuery * QUERY_SIZE_IN_WORD * sizeof(int),cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemcpy(d_queryIdBuf_temp,d_queryIdBuf,
				numQuery*sizeof(int), cudaMemcpyDeviceToDevice));

		cutilSafeCall(cudaMemset(d_queryBuf, 0, numQuery*QUERY_SIZE_IN_WORD*sizeof(int)));
		cutilSafeCall(cudaMemset(d_queryIdBuf, 0, numQuery*sizeof(int)));
                d_queryBufExec = d_queryBuf_temp;
                d_queryIdBufExec = d_queryIdBuf_temp;
	}

	numBlock = 512; 
	numThread = 128; 
	passId = 0;
	numUndone = numQuery;
#if SORT_QUERY
	d_queryBufExec = d_Queries_Fwd;
#else
	d_queryBufExec = d_Queries;
#endif
	d_queryIdBufExec = d_QueryId;
	while(numUndone > 0) {
		if(numUndone < SMALL_DATA_SIZE){
			numBlock = SMALL_BLOCK_SIZE;
			numThread = SMALL_THREAD_SIZE;
		}
		maxNumResultPerThread = getNumOutputPerThread(gpuResultBufSize, numBlock, numThread, numUndone);
		cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
		cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
		passId++;
		sharedMemSize = numThread * OCC_COUNT_PER_THREAD * sizeof(int);
		BWT2MismatchSearchCaseC_shareBuf_kernel<<<numBlock, numThread, sharedMemSize>>>(d_result, 
				d_numResultPerQuery,maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_bwt, d_revbwt); 
		cutilCheckMsg("BWT1Mismatch forward_kernel 1"); 
		
		uint64 numAlignment1MisForwardPass1 = 0; 
		int numUndone1MisForwardPass1 = 0; 

		collectUndoneGPU_shareBuf(       &numAlignment1MisForwardPass1, 
				d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
				d_result, d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_temp, d_queryLabel, numBlock,numThread, h_result, h_resultNum); 
		
		printf("pass %d caseC undone num = %lu\n",passId, numUndone1MisForwardPass1); 
		printf("numAlignmentCaseC= %lu\n\n", numAlignment1MisForwardPass1); 

		numAlignment += numAlignment1MisForwardPass1; 
		numUndone = numUndone1MisForwardPass1;

		cutilSafeCall(cudaMemcpy(d_queryBuf_temp,d_queryBuf, 
				numQuery * QUERY_SIZE_IN_WORD * sizeof(int),cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemcpy(d_queryIdBuf_temp,d_queryIdBuf,
				numQuery*sizeof(int), cudaMemcpyDeviceToDevice));

		cutilSafeCall(cudaMemset(d_queryBuf, 0, numQuery*QUERY_SIZE_IN_WORD*sizeof(int)));
		cutilSafeCall(cudaMemset(d_queryIdBuf, 0, numQuery*sizeof(int)));
                d_queryBufExec = d_queryBuf_temp;
		d_queryIdBufExec = d_queryIdBuf_temp;
	}

	numBlock = 512; 
	numThread = 128; 
	passId = 0;
	numUndone = numQuery;
#if SORT_QUERY
	d_queryBufExec = d_Queries_Fwd;
#else
	d_queryBufExec = d_Queries;
#endif
	d_queryIdBufExec = d_QueryId;
	while(numUndone > 0) {
		if(numUndone < SMALL_DATA_SIZE){
			numBlock = SMALL_BLOCK_SIZE;
			numThread = SMALL_THREAD_SIZE;
		}
		maxNumResultPerThread = getNumOutputPerThread(gpuResultBufSize, numBlock, numThread, numUndone);
		cutilSafeCall(cudaMemset(d_numResultPerQuery, 0 ,numQuery*sizeof(int2)));
		cutilSafeCall(cudaMemset(d_result, 0 , gpuResultBufSize));
		passId++;
		sharedMemSize = numThread * OCC_COUNT_PER_THREAD * sizeof(int);
		BWT2MismatchSearchCaseD_shareBuf_kernel<<<numBlock, numThread,sharedMemSize>>>(d_result, 
				d_numResultPerQuery,maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_bwt, d_revbwt); 
		cutilCheckMsg("BWT1Mismatch forward_kernel 1"); 
		
		uint64 numAlignment1MisForwardPass1 = 0; 
		int numUndone1MisForwardPass1 = 0; 

		collectUndoneGPU_shareBuf(       &numAlignment1MisForwardPass1, 
				d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
				d_result, d_numResultPerQuery, maxNumResultPerThread, 
				d_queryBufExec, d_queryIdBufExec, numUndone, d_temp, d_queryLabel, numBlock,numThread, h_result, h_resultNum); 
		
		printf("pass %d caseD undone num = %lu\n",passId, numUndone1MisForwardPass1); 
		printf("numAlignmentCaseC= %lu\n\n", numAlignment1MisForwardPass1); 

		numAlignment += numAlignment1MisForwardPass1; 
		numUndone = numUndone1MisForwardPass1;
		cutilSafeCall(cudaMemcpy(d_queryBuf_temp,d_queryBuf, 
				numQuery * QUERY_SIZE_IN_WORD * sizeof(int),cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemcpy(d_queryIdBuf_temp,d_queryIdBuf,
				numQuery*sizeof(int), cudaMemcpyDeviceToDevice));
		cutilSafeCall(cudaMemset(d_queryBuf, 0, numQuery*QUERY_SIZE_IN_WORD*sizeof(int)));
		cutilSafeCall(cudaMemset(d_queryIdBuf, 0, numQuery*sizeof(int)));
                d_queryBufExec = d_queryBuf_temp;
		d_queryIdBufExec = d_queryIdBuf_temp;
	}




	const int numMatched = getNumMatchedQuery(d_queryLabel, queryNum); 
	printf("===> numAlignment: %lu\n", numAlignment); 
	printf("===> #matched query: %d\n", numMatched); 

#if SORT_QUERY 
	free(sortedQueryForward); 
	free(sortedQueryBackward); 
#endif 
	return numAlignment; 



}

uint64 GPU2BWTSearch_noConcurent(const unsigned int * h_Queries, unsigned int queryNum,  
		BWT *d_bwt, BWT *d_revbwt, 
		st_MemoryBuf* GPUMemBuf) 
{ 

	PRINT_FUNC_NAME; 

	/*GPU memory*/ 
	unsigned int numQuery = queryNum; 
	const uint64 gpuResultBufSize = 2200000000; //just used for result output per kernel 
	uint3* d_result = (uint3*)GPUMemBuf->newMem(gpuResultBufSize); 
	int* d_numResultPerQuery = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery); 
	unsigned int* d_queryBuf = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
	int* d_queryIdBuf = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery); 
	int* d_temp = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery); 


	/***********************Sort query***********************/ 
#if SORT_QUERY   
	unsigned int * sortedQueryForward = SortQueriesForward(h_Queries, numQuery); 
	unsigned int * sortedQueryBackward = SortQueriesBackward(h_Queries, numQuery); 
	unsigned int* d_Queries_Fwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
	unsigned int* d_Queries_Bwd = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
	TOGPU(d_Queries_Fwd, sortedQueryForward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
	TOGPU(d_Queries_Bwd, sortedQueryBackward, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
#else    
	unsigned int* d_Queries = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
	TOGPU(d_Queries, h_Queries, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery); 
#endif 
	int* origQueryId = (int*)malloc(sizeof(int)*queryNum); 
	for(int i = 0; i < queryNum; i++) { 
		origQueryId[i] = i; 
	} 
	int* d_QueryId = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum); 
	TOGPU(d_QueryId, origQueryId, sizeof(int)*queryNum); 
	int* d_queryLabel = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum); 
	cutilSafeCall(cudaMemset(d_queryLabel, 0, sizeof(int)*queryNum)); 
	printf("GPU memory done.\n"); 



	/*kernel configuration and others*/ 
#ifdef PREFER_L1 
	cudaFuncSetCacheConfig(exactMatch_kernel, cudaFuncCachePreferL1); 
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferL1); 
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferL1); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferL1); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferL1); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferL1); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferL1); 
#else  
	cudaFuncSetCacheConfig(exactMatch_kernel, cudaFuncCachePreferShared); 
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferShared); 
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferShared); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferShared); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferShared); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferShared); 
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferShared); 

#endif 
	int maxNumResultPerQuery; 
	int numBlock, numThread; 
	int sharedMemSize = 0; 
	uint64 numAlignment = 0; 
	sharedMemSize = 0; //sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numThread; 


	/*******************************************************/ 
	/*exact match pass1 only*/ 
	printTitle("exact match"); 
	numBlock = 512; 
	//numBlock = 1204; 
	numThread = 128; 
#if SORT_QUERY 
	exactMatch_kernel<<<numBlock, numThread>>>(d_result,  
			d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt); 
#else 
	exactMatch_kernel<<<numBlock, numThread>>>(d_result,  
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt); 
#endif 
	cutilCheckMsg("BWT exact match kernel 1"); 
	uint3* resultExact = (uint3 *)malloc(sizeof(uint3) * numQuery ); 
	cutilSafeCall(cudaMemcpy(resultExact, d_result, sizeof(uint3)*numQuery, cudaMemcpyDeviceToHost)); 

	for(int i=0; i<numQuery; i++) { 
		if(resultExact[i].x  < numQuery) 
			numAlignment += resultExact[i].z - resultExact[i].y + 1;         
	} 
	printf("numAlignment exact match: %lu\n\n", numAlignment); 


	/*******************************************************/ 
	/*1 mismatch*/ 
	/*1 mismatch forward*/ 
	//pass 1 
	printTitle("1 mismatch"); 
	printf("1 mismatch forward search: %d\n", numQuery); 
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery); 
	numBlock = 512; 
	//      numBlock = 1204; 
	numThread = 128; 
#if SORT_QUERY 
	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery, 
			d_Queries_Fwd, numQuery, d_bwt, d_revbwt); 
#else 
	BWT1MismatchForward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery, 
			d_Queries, numQuery, d_bwt, d_revbwt); 
#endif   
	cutilCheckMsg("BWT1Mismatch forward_kernel 1"); 
	uint64 numAlignment1MisForwardPass1 = 0; 
	int numUndone1MisForwardPass1 = 0; 

#if SORT_QUERY 
	collectUndoneGPU(       &numAlignment1MisForwardPass1, d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel); 
#else    
	collectUndoneGPU(       &numAlignment1MisForwardPass1, d_queryBuf, d_queryIdBuf, &numUndone1MisForwardPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel); 
#endif   
	printf("1MisForward undone num = %lu\n", numUndone1MisForwardPass1); 
	printf("numAlignment1MisForwardPass1 = %lu\n\n", numAlignment1MisForwardPass1); 

	numAlignment += numAlignment1MisForwardPass1; 
	//pass 2 
	if(numUndone1MisForwardPass1 > 0) {
		maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisForwardPass1); 

		BWT1MismatchForward_Kernel<<<numBlock, numThread>>> 
			(d_result, d_numResultPerQuery, 
			 maxNumResultPerQuery, d_queryBuf, numUndone1MisForwardPass1, d_bwt, d_revbwt); 

		uint64 numAlignment1MisForwardPass2 = 0; 
		int numUndone1MisForwardPass2 = 0; 
		countAlignGPU(  &numAlignment1MisForwardPass2, d_result, d_numResultPerQuery, 
				maxNumResultPerQuery, numUndone1MisForwardPass1, 
				d_temp, d_queryIdBuf, d_queryLabel); 
		printf("numUndone1MisForwardPass2 = %d\n", numUndone1MisForwardPass2); 
		assert(numUndone1MisForwardPass2 == 0); 
		printf("numAlignment1MisForwardPass2 = %lu\n\n", numAlignment1MisForwardPass2); 
		numAlignment += numAlignment1MisForwardPass2; 
	}




	/* backward */ 
	//pass 1 
	printf("1 mismatch backward search: %d\n", numQuery); 
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery); 
	numBlock = 512; 
	//        numBlock = 1204; 
	numThread = 128; 
#if SORT_QUERY 
	BWT1MismatchBackward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery, 
			d_Queries_Bwd, numQuery, d_bwt, d_revbwt); 
#else 
	BWT1MismatchBackward_Kernel<<<numBlock, numThread>>>(d_result, d_numResultPerQuery,maxNumResultPerQuery, 
			d_Queries, numQuery, d_bwt, d_revbwt); 
#endif 
	cutilCheckMsg("BWT1Mismatch backward_kernel 1"); 


	uint64 numAlignment1MisBackwardPass1 = 0; 
	int numUndone1MisBackwardPass1 = 0; 
#if SORT_QUERY 
	collectUndoneGPU(       &numAlignment1MisBackwardPass1, d_queryBuf, d_queryIdBuf, &numUndone1MisBackwardPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries_Bwd, d_QueryId, queryNum, d_temp, d_queryLabel); 
#else 
	collectUndoneGPU(       &numAlignment1MisBackwardPass1, d_queryBuf, d_queryIdBuf, &numUndone1MisBackwardPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel); 
#endif   
	printf("1MisBackward undone num = %lu\n", numUndone1MisBackwardPass1); 
	printf("numAlignment1MisBackwardPass1 = %lu\n\n", numAlignment1MisBackwardPass1); 
	numAlignment += numAlignment1MisBackwardPass1; 
	//pass 2 
	if(numUndone1MisBackwardPass1 > 0) {
		const int maxNumResultPerQuery1B = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisBackwardPass1); 
		BWT1MismatchBackward_Kernel<<<numBlock, numThread, 0>>> 
			(d_result, d_numResultPerQuery, 
			 maxNumResultPerQuery, d_queryBuf, numUndone1MisBackwardPass1, d_bwt, d_revbwt); 
		uint64 numAlignment1MisBackwardPass2 = 0; 
		int numUndone1MisBackwardPass2 = 0; 
		countAlignGPU(  &numAlignment1MisBackwardPass2, d_result, d_numResultPerQuery, 
				maxNumResultPerQuery, numUndone1MisBackwardPass1, 
				d_temp, d_queryIdBuf, d_queryLabel); 
		printf("numUndone1MisBackwardPass2 = %d\n", numUndone1MisBackwardPass2); 
		assert(numUndone1MisBackwardPass2 == 0); 
		printf("numAlignment1MisBackwardPass2 = %lu\n\n", numAlignment1MisBackwardPass2); 
		numAlignment += numAlignment1MisBackwardPass2; 
	} 




	/*******************************************************/ 
	/*two mismatches*/ 
	printTitle("2 mismatches"); 
	//case A, pass 1 
	printf("Case A, pass 1, numQuery: %d\n", numQuery); 
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery); 
#if SORT_QUERY 
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize>>> 
		(d_result, d_numResultPerQuery,  
		 maxNumResultPerQuery, d_Queries_Bwd, d_QueryId, numQuery, d_bwt, d_revbwt); 
#else 
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize>>> 
		(d_result, d_numResultPerQuery,  
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseA_kernel 1");
	uint64 numAlignmentCaseAPass1 = 0;
	int numUndoneCaseAPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignmentCaseAPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseAPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Bwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignmentCaseAPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseAPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseAPass1 = %lu\n", numAlignmentCaseAPass1);
	printf("numUndoneCaseAPass1 = %lu\n\n", numUndoneCaseAPass1);
	numAlignment += numAlignmentCaseAPass1; 

	numBlock = 448;
	numThread = 64;
	const int maxNumResultPerQueryA2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseAPass1);
	printf("pass 2: numBlock = %d, numThread = %d, maxNumPerQuery = %d\n", numBlock, numThread, maxNumResultPerQueryA2);
	BWT2MismatchSearchCaseA_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQueryA2, d_queryBuf, d_queryIdBuf,
		 numUndoneCaseAPass1, d_bwt, d_revbwt);
	printf("Case A, pass 2, #query = %d\n", numUndoneCaseAPass1);
	uint64 numAlignmentCaseAPass2 = 0;
	countAlignGPU(  &numAlignmentCaseAPass2, d_result, d_numResultPerQuery,
			maxNumResultPerQueryA2, numUndoneCaseAPass1,
			d_temp, d_queryIdBuf, d_queryLabel);
	printf("numAlignmentCaseAPass2 = %lu\n\n", numAlignmentCaseAPass2);
	numAlignment += numAlignmentCaseAPass2;		



	//case B, pass 1

	numBlock = 512; 
	numThread = 128;
	numQuery = queryNum;
	printf("Case B, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#if SORT_QUERY
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries_Fwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseB_kernel 1");
	uint64 numAlignmentCaseBPass1 = 0;
	int numUndoneCaseBPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignmentCaseBPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseBPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignmentCaseBPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseBPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseBPass1 = %lu\n", numAlignmentCaseBPass1);
	printf("numUndoneCaseBPass1 = %lu\n\n", numUndoneCaseBPass1);
	numAlignment += numAlignmentCaseBPass1; 

	numBlock = 448;
	numThread = 64;
	const int maxNumResultPerQueryB2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseBPass1);
	printf("pass 2: numBlock = %d, numThread = %d, maxNumPerQuery = %d\n", numBlock, numThread, maxNumResultPerQueryB2);
	BWT2MismatchSearchCaseB_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQueryB2, d_queryBuf, d_queryIdBuf,
		 numUndoneCaseBPass1, d_bwt, d_revbwt);
	printf("Case B, pass 2, #query = %d\n", numUndoneCaseBPass1);
	uint64 numAlignmentCaseBPass2 = 0;
	countAlignGPU(  &numAlignmentCaseBPass2, d_result, d_numResultPerQuery,
			maxNumResultPerQueryB2, numUndoneCaseBPass1,
			d_temp, d_queryIdBuf, d_queryLabel);
	printf("numAlignmentCaseBPass2 = %lu\n\n", numAlignmentCaseBPass2);
	numAlignment += numAlignmentCaseBPass2;







	//case C, pass 1
	numBlock = 512; 
	numThread = 128;
	numQuery = queryNum;
	printf("Case C, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#if SORT_QUERY
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries_Fwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseC_kernel 1");
	uint64 numAlignmentCaseCPass1 = 0;
	int numUndoneCaseCPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignmentCaseCPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseCPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignmentCaseCPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseCPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseCPass1 = %lu\n", numAlignmentCaseCPass1);
	printf("numUndoneCaseCPass1 = %lu\n\n", numUndoneCaseCPass1);
	numAlignment += numAlignmentCaseCPass1;

	numBlock = 448;
	numThread = 64;
	const int maxNumResultPerQueryC2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseCPass1);
	printf("pass 2: numBlock = %d, numThread = %d, maxNumPerQuery = %d\n", numBlock, numThread, maxNumResultPerQueryC2);
	BWT2MismatchSearchCaseC_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQueryC2, d_queryBuf, d_queryIdBuf,
		 numUndoneCaseCPass1, d_bwt, d_revbwt);
	printf("Case C, pass 2, #query = %d\n", numUndoneCaseCPass1);
	uint64 numAlignmentCaseCPass2 = 0;
	countAlignGPU(  &numAlignmentCaseCPass2, d_result, d_numResultPerQuery,
			maxNumResultPerQueryC2, numUndoneCaseCPass1,
			d_temp, d_queryIdBuf, d_queryLabel);
	printf("numAlignmentCaseCPass2 = %lu\n\n", numAlignmentCaseCPass2);
	numAlignment += numAlignmentCaseCPass2;

	//case D, pass 1
	numBlock = 512; 
	numThread = 128;
	numQuery = queryNum;
	printf("Case D, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
#if SORT_QUERY
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries_Fwd, d_QueryId, numQuery, d_bwt, d_revbwt);
#else
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
#endif
	cutilCheckMsg("BWT2MismatchSearchCaseD_kernel 1");
	uint64 numAlignmentCaseDPass1 = 0;
	int numUndoneCaseDPass1 = 0;
#if SORT_QUERY
	collectUndoneGPU(       &numAlignmentCaseDPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseDPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries_Fwd, d_QueryId, queryNum, d_temp, d_queryLabel);
#else
	collectUndoneGPU(       &numAlignmentCaseDPass1, d_queryBuf, d_queryIdBuf, &numUndoneCaseDPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_temp, d_queryLabel);
#endif
	printf("numAlignmentCaseDPass1 = %lu\n", numAlignmentCaseDPass1);
	printf("numUndoneCaseDPass1 = %lu\n\n", numUndoneCaseDPass1);
	numAlignment += numAlignmentCaseDPass1;

	numBlock = 448;
	numThread = 64;
	const int maxNumResultPerQueryD2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseDPass1);
	printf("pass 2: numBlock = %d, numThread = %d, maxNumPerQuery = %d\n", numBlock, numThread, maxNumResultPerQueryD2);
	BWT2MismatchSearchCaseD_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_result, d_numResultPerQuery,
		 maxNumResultPerQueryD2, d_queryBuf, d_queryIdBuf,
		 numUndoneCaseDPass1, d_bwt, d_revbwt);
	printf("Case D, pass 2, #query = %d\n", numUndoneCaseDPass1);
	uint64 numAlignmentCaseDPass2 = 0;
	countAlignGPU(  &numAlignmentCaseDPass2, d_result, d_numResultPerQuery,
			maxNumResultPerQueryD2, numUndoneCaseDPass1,
			d_temp, d_queryIdBuf, d_queryLabel);
	printf("numAlignmentCaseDPass2 = %lu\n\n", numAlignmentCaseDPass2);
	numAlignment += numAlignmentCaseDPass2;



	const int numMatched = getNumMatchedQuery(d_queryLabel, queryNum); 
	printf("===> numAlignment: %lu\n", numAlignment); 
	printf("===> #matched query: %d\n", numMatched); 

#if SORT_QUERY 
	free(sortedQueryForward); 
	free(sortedQueryBackward); 
#endif 
	return numAlignment; 
}

uint64 GPU2BWTSearchTwoStrandLimit(long* d_numAlign,
		const unsigned int* h_Queries, const unsigned int* h_revQueries,
		unsigned int queryNum,
		BWT *d_bwt, BWT *d_revbwt,
		st_MemoryBuf* GPUMemBuf, uint3* h_result, uint32* h_resultNum) {

	PRINT_FUNC_NAME;

	cutilSafeCall(cudaMemset(d_numAlign, 0, sizeof(long)*queryNum));
	GPUMemBuf->reset();

	//the positive strand
	GPU2BWTSearchLimit(d_numAlign, h_Queries, queryNum, d_bwt, d_revbwt, GPUMemBuf, h_result, h_resultNum);
	GPUMemBuf->reset();

	//the negative strand
	GPU2BWTSearchLimit(d_numAlign, h_revQueries, queryNum, d_bwt, d_revbwt, GPUMemBuf, h_result, h_resultNum);
	GPUMemBuf->reset();

	return 0;
}

uint64 GPU2BWTSearchLimit(long* d_numAlign,
		const unsigned int * h_Queries, unsigned int queryNum, 
		BWT *d_bwt, BWT *d_revbwt,
		st_MemoryBuf* GPUMemBuf, uint3* h_result, uint32* h_resultNum) {

	PRINT_FUNC_NAME;


	/*GPU memory*/
	int numQuery = queryNum;
	const uint64 gpuResultBufSize = 450000000; //just used for result output per kernel
	uint3* d_result = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	int* d_numResultPerQuery = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	///now for the concurrent kernels
	unsigned int* d_queryBufA = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBufB = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBufC = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBufD = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBuf1F = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	unsigned int* d_queryBuf1B = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* d_queryIdBufA = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBufB = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBufC = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBufD = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBuf1F = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_queryIdBuf1B = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
/*	uint3* d_resultA = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_resultB = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_resultC = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_resultD = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_result1F = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
	uint3* d_result1B = (uint3*)GPUMemBuf->newMem(gpuResultBufSize);
*/	int* d_numResultPerQueryA = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQueryB = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQueryC = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQueryD = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQuery1F = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_numResultPerQuery1B = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	int* d_temp = (int*)GPUMemBuf->newMem(sizeof(int)*numQuery);
	unsigned int* d_Queries = (unsigned int*)GPUMemBuf->newMem(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	TOGPU(d_Queries, h_Queries, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	int* origQueryId = (int*)malloc(sizeof(int)*queryNum);
	for(int i = 0; i < queryNum; i++) {
		origQueryId[i] = i;
	}
	int* d_QueryId = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	TOGPU(d_QueryId, origQueryId, sizeof(int)*queryNum);
	int* d_queryLabel = (int*)GPUMemBuf->newMem(sizeof(int)*queryNum);
	cutilSafeCall(cudaMemset(d_queryLabel, 0, sizeof(int)*queryNum));
	printf("GPU memory done.\n");



	/*kernel configuration and others*/
	cudaFuncSetCacheConfig(BWT1MismatchForward_Kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT1MismatchBackward_Kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseA_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseB_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseC_kernel, cudaFuncCachePreferL1);
	cudaFuncSetCacheConfig(BWT2MismatchSearchCaseD_kernel, cudaFuncCachePreferL1);
	int maxNumResultPerQuery;
	int numBlock, numThread;
	int sharedMemSize = 0;
	uint64 numAlignment = 0;
	int nstreams = 6;
	cudaStream_t *streams = (cudaStream_t*) malloc(nstreams * sizeof(cudaStream_t));
	for(int i = 0; i < nstreams; i++) {
		cutilSafeCall( cudaStreamCreate(&(streams[i])) );
	}
	numBlock = 128;
	numThread = 64;
	sharedMemSize = 0; //sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numThread;

	/*******************************************************/
	/*exact match pass1 only*/

	printTitle("exact match");
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	exactMatch_limit_kernel<<<numBlock, numThread>>>(d_numAlign, d_result,d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
	cutilCheckMsg("BWT exact match kernel 1");
	uint3* resultExact = (uint3 *)malloc(sizeof(uint3) * numQuery);
	cutilSafeCall(cudaMemcpy(resultExact, d_result, sizeof(uint3)*numQuery, cudaMemcpyDeviceToHost));

	printf("here\n");
	for(int i=0; i<numQuery; i++) {
		if(resultExact[i].x  < numQuery){
			numAlignment += resultExact[i].z - resultExact[i].y + 1;
			h_result[*h_resultNum] = resultExact[i];
			(*h_resultNum)++;
		}
	}
	printf("numAlignment exact match: %lu\n\n", numAlignment);



	/**********************************************************/
	/*1 mismatch*/

	printTitle("1 mismatch");
	printf("1 mismatch forward search: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	uint64 numAlignment1MisForwardPass1 = 0;
	unsigned int timer = 0;
	int numUndone1MisForwardPass1 = 0;
	uint64 numAlignment1MisBackwardPass1 = 0;
	int numUndone1MisBackwardPass1 = 0;
	startTimer(&timer);
	BWT1MismatchForward_limit_Kernel<<<numBlock, numThread>>>(d_numAlign, d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);

	cutilCheckMsg("BWT1Mismatch forward_kernel 1");
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer, "1 mismatch, forward, kernel");
	collectUndoneGPULimit(d_numAlign, d_queryBuf1F, d_queryIdBuf1F, &numUndone1MisForwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_queryLabel, d_temp, h_result, h_resultNum);
	printf("1MisForward undone num = %lu\n", numUndone1MisForwardPass1);
	printf("numAlignment1MisForwardPass1 = %lu\n\n", numAlignment1MisForwardPass1);

	numAlignment += numAlignment1MisForwardPass1;

	printf("1 mismatch backward search: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	timer = 0;
	startTimer(&timer);
	BWT1MismatchBackward_limit_Kernel<<<numBlock, numThread>>>(d_numAlign, d_result, d_numResultPerQuery,maxNumResultPerQuery,
			d_Queries, d_QueryId, numQuery, d_bwt, d_revbwt);
	cutilCheckMsg("BWT1Mismatch backward_kernel 1");
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer, "1 mismatch, backward, kernel");
	collectUndoneGPULimit(d_numAlign, d_queryBuf1B, d_queryIdBuf1B, &numUndone1MisBackwardPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery,
			d_Queries, d_QueryId, queryNum, d_queryLabel, d_temp, h_result, h_resultNum);
	printf("1MisBackward undone num = %lu\n", numUndone1MisBackwardPass1);
	printf("numAlignment1MisBackwardPass1 = %lu\n\n", numAlignment1MisBackwardPass1);
	numAlignment += numAlignment1MisBackwardPass1;

	//////////////////////


	/*******************************************************/
	printTitle("2 mismatches");
	//pass 1

	//case A, pass 1
	printf("Case A, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	BWT2MismatchSearchCaseALimit_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_numAlign, d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, 
		 d_bwt, d_revbwt);
	cutilCheckMsg("BWT2MismatchSearchCaseA_kernel 1");
	int numUndoneCaseAPass1 = 0;
	collectUndoneGPULimit(d_numAlign, d_queryBufA, d_queryIdBufA, &numUndoneCaseAPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_queryLabel, d_temp, h_result, h_resultNum);
	printf("case A done.\n");



	//case B, pass 1
	numQuery = queryNum;
	printf("Case B, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	BWT2MismatchSearchCaseBLimit_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_numAlign, d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, 
		 d_bwt, d_revbwt);
	cutilCheckMsg("BWT2MismatchSearchCaseB_kernel 1");
	int numUndoneCaseBPass1 = 0;
	collectUndoneGPULimit(d_numAlign, d_queryBufB, d_queryIdBufB, &numUndoneCaseBPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_queryLabel, d_temp, h_result, h_resultNum);
	printf("case B done.\n");


	//case C, pass 1
	numQuery = queryNum;
	printf("Case C, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	BWT2MismatchSearchCaseCLimit_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_numAlign, d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, 
		 d_bwt, d_revbwt);
	cutilCheckMsg("BWT2MismatchSearchCaseC_kernel 1");
	int numUndoneCaseCPass1 = 0;
	collectUndoneGPULimit(d_numAlign, d_queryBufC, d_queryIdBufC, &numUndoneCaseCPass1, 
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_queryLabel, d_temp, h_result, h_resultNum);

	printf("case C done.\n");


	//case D, pass 1
	numQuery = queryNum;
	printf("Case D, pass 1, numQuery: %d\n", numQuery);
	maxNumResultPerQuery = getNumOutputPerQuery(gpuResultBufSize, numQuery);
	BWT2MismatchSearchCaseDLimit_kernel<<<numBlock, numThread, sharedMemSize>>>
		(d_numAlign, d_result, d_numResultPerQuery, 
		 maxNumResultPerQuery, d_Queries, d_QueryId, numQuery, 
		 d_bwt, d_revbwt);
	cutilCheckMsg("BWT2MismatchSearchCaseD_kernel 1");
	int numUndoneCaseDPass1 = 0;
	collectUndoneGPULimit(d_numAlign, d_queryBufD, d_queryIdBufD, &numUndoneCaseDPass1,
			d_result, d_numResultPerQuery, maxNumResultPerQuery, 
			d_Queries, d_QueryId, queryNum, d_queryLabel, d_temp, h_result, h_resultNum);
	printf("case D done.\n");




	/////////////////////////////////////////////////////
	//Pass 2
/*
	numBlock = 448;
	numThread = 64;
	sharedMemSize = 0; //sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numThread;

	/////////////////////
	//pass 2, stream execution
	printf("pass 2: numBlock = %d, numThread = %d\n", numBlock, numThread);
	const int maxNumResultPerQuery1F = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisForwardPass1);
	const int maxNumResultPerQuery1B = getNumOutputPerQuery(gpuResultBufSize, numUndone1MisBackwardPass1);
	const int maxNumResultPerQueryA2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseAPass1);
	const int maxNumResultPerQueryB2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseBPass1);
	const int maxNumResultPerQueryC2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseCPass1);
	const int maxNumResultPerQueryD2 = getNumOutputPerQuery(gpuResultBufSize, numUndoneCaseDPass1);




	//kernels
	//1 mismatch

	BWT1MismatchForward_limit_Kernel<<<numBlock, numThread, 0, streams[0]>>>
		(d_numAlign,d_result1F, d_numResultPerQuery1F,
		 maxNumResultPerQuery1F, d_queryBuf1F,d_queryIdBuf1F, numUndone1MisForwardPass1, d_bwt, d_revbwt);

	BWT1MismatchBackward_limit_Kernel<<<numBlock, numThread, 0, streams[0]>>>
		(d_numAlign, d_result1B, d_numResultPerQuery1B,
		 maxNumResultPerQuery1B, d_queryBuf1B, d_queryIdBuf1B, numUndone1MisBackwardPass1, d_bwt, d_revbwt);



	//two mismatches
	BWT2MismatchSearchCaseALimit_kernel<<<numBlock, numThread, sharedMemSize, streams[0]>>>
		(d_numAlign, d_resultA, d_numResultPerQueryA, 
		 maxNumResultPerQueryA2, d_queryBufA, d_queryIdBufA, 
		 numUndoneCaseAPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseBLimit_kernel<<<numBlock, numThread, sharedMemSize, streams[0]>>>
		(d_numAlign, d_resultB, d_numResultPerQueryB, 
		 maxNumResultPerQueryB2, d_queryBufB, d_queryIdBufB, 
		 numUndoneCaseBPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseCLimit_kernel<<<numBlock, numThread, sharedMemSize, streams[0]>>>
		(d_numAlign, d_resultC, d_numResultPerQueryC, 
		 maxNumResultPerQueryC2, d_queryBufC, d_queryIdBufC, 
		 numUndoneCaseCPass1, d_bwt, d_revbwt);
	BWT2MismatchSearchCaseDLimit_kernel<<<numBlock, numThread, sharedMemSize, streams[0]>>>
		(d_numAlign, d_resultD, d_numResultPerQueryD, 
		 maxNumResultPerQueryD2, d_queryBufD, d_queryIdBufD, 
		 numUndoneCaseDPass1, d_bwt, d_revbwt);
*/
	return numAlignment;
}





/**
  @return the number of alignment
 */


#endif
